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SYNTHETIC APERTURE IMAGING OF DIRECTION AND 
FREQUENCY DEPENDENT REFLECTIVITIES 


LILIANA BORCEA*, MIGUEL MOSCOSO t, GEORGE PAPANICOLAOUt, AND 

GHRYSOULA TSOGKAS 

Abstract. We introduce a synthetic aperture imaging framework that takes into consideration 
directional dependence of the reflectivity that is to be imaged, as well as its frequency dependence. 
We use an minimization approach that is coordinated with data segmentation so as to fuse 
information from multiple sub-apertures and frequency sub-bands. We analyze this approach from 
first principles and assess its performance with numerical simulations in an X-band radar regime. 

Key words, synthetic aperture imaging, reflectivity, minimal support optimization. 

1. Introduction. We introduce and analyze a novel algorithm for synthetic 
aperture radar (SAR) imaging, where a moving receive-transmit platform probes a 
remote region with signals f(t) and records the scattered waves. The platform spans a 
large synthetic aperture so that high resolution images of the region may be obtained 
by processing the recorded data. A related application is inverse synthetic aperture 
radar (ISAR), where the receive-transmit antenna is stationary, and the synthetic 
aperture is due to the motion of an unknown scatterer. If this motion is known or 
can be estimated, the problem can be restated mathematically as SAR imaging of the 
scatterer, using the reference frame that moves with it. 

A schematic of the SAR imaging setup is in Figure 1.1. The recordings u{s^t) at 
the moving receive-transmit platform depend on two time variables: the slow time s 
and the fast time t. The slow time parametrizes the trajectory of the platform, and 
it is discretized in uniform steps called the pulse repetition rate. At time s the 
platform is at location r{s). It emits the signal f{t) and receives the backscattered 
returns u{s,t). The fast time t runs between consecutive signal emissions t G (0,hs), 
and we assume a separation of time scales: The duration of f{t) is smaller than the 
round trip travel time of the waves between the sensor and the imaging region, and 
the latter is smaller than hg. 

In the usual synthetic aperture image formulation the reflectivity is modeled as a 
two dimensional function of location y on a surface of known topography, say flat for 
simplicity. The assumption is that each point on the surface reflects the waves the 
same way in all directions, independent of the direction and frequency of the incident 
waves. This simplifies the imaging process and makes the inverse problem formally 
determined: the data are two-dimensional and so is the unknown reflectivity function. 

The reflectivity can be reconstructed by the reverse time migration formula [20, 
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Fig. 1.1. Setup for imaging with a synthetic aperture. 


Here Sj are the slow time emission-recording instants, spaced by and the im¬ 
age is formed by superposing over the platform trajectory the data u{sj,t)^ match- 
filtered with the time reversed emitted signal /(t), delayed by the roundtrip travel 
time 2r(sj,y) = 2\T{sj) — y|/c between the platform location T{sj) and the imaging 
point y. The bar denotes complex conjugate and c is the wave speed in the medium 
which is assumed homogeneous. 

The assumption of an isotropic reflectivity may not always be justified in applica¬ 
tions. Backscatter reflectivities are in general functions of five variables: the location 
y on the known (flat) surface, the two angles of incidence and the frequency. Thus, 
the inverse problem is underdetermined and we cannot expect a reconstruction of the 
five dimensional reflectivity with a migration approach. Direct application of (1.1) 
will produce low-resolution images of some effective, position-dependent reflectivity, 
and there will be no information about the directivity and frequency dependence of 
the actual reflectivity. 

The reconstruction of frequency dependent reflectivities with synthetic aperture 
radar has been considered in [8], where Doppler effects are shown to be useful in 
inversion, and in [21, 11], where data are segmented over frequency sub-bands, and 
then images are formed separately, for each data subset. Data segmentation is a 
natural idea, and we show here how to use it for reconstructing both frequency and 
direction dependent reflectivities. 

The main result in this paper is the introduction and analysis of an algorithm for 
imaging direction and frequency dependent reflectivities of strong, localized scatterers. 
This algorithm is based on optimization. It reconstructs reflectivities of localized 
scatterers by seeking among all those that fit the data model the ones with minimal 
spatial support. Array imaging algorithms based on £i optimization are proposed and 
analyzed in [1, 18, 5, 4, 14, 13, 2]. They consider only isotropic, frequency independent 
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reflectivities. 

A direct extension of optimization methods to imaging direction and frequency 
dependent reflectivities amounts to solving a grand optimization problem for a very 
long vector p of unknowns, the discretized reflectivity over spatial locations on the 
imaging grid, the angles of incidence/backscatter and the frequency. It has consider¬ 
able computational complexity because of the high dimension of the space in which 
the discretized reflectivity vector lies. It also does not take into account the fact that 
many unknowns are tied to the same spatial location points within the discretized 
image window. 

The synthetic aperture imaging algorithm introduced in this paper is designed to 
reconstruct efficiently direction and frequency dependent reflectivities by combining 
two main ideas: The first is to divide the data over carefully calibrated sub-apertures 
and frequency sub-bands, and solve an optimization problem to estimate the reflec¬ 
tivity for each data subset. Data segmentation is useful assuming that the reflectivity 
changes continuously with the direction of probing and the frequency, so that we 
can approximate it by a piecewise constant function, pointwise in the imaging win¬ 
dow. Over a sub-aperture of small enough linear size a, the platform receives scattered 
waves from a narrow cone with opening angle of the order a/L, where L is the distance 
from the platform to the imaging window, and we can approximate the reflectivity 
by that at the center angle. Similarly, we can approximate the reflectivity by a con¬ 
stant over a small enough frequency sub-band. Then, we can use optimization 
to estimate the reflectivity as a function of location for each data subset. The size 
of the sub-apertures and sub-bands determine the resolution of the reconstruction. 
The larger they are, the better the expected spatial resolution of the reflectivity. But 
the resolution is worse over direction and frequency dependence. The calibration of 
the data segmentation over sub-apertures and sub-bands reflects this trade-off. The 
second idea combines the optimizations by seeking reflectivities that have common 
spatial support. Instead of a single vector p, the unknown is a matrix with columns 
of spatially discretized reflectivities. Each column corresponds to a direction of prob¬ 
ing from a sub-aperture and a central frequency in a sub-band. The values of the 
entries in the columns are different, but they are zero (negligible) in the same rows. 
Moreover, the forward model, which is derived here from first principles, maps each 
column of the reflectivity matrix to the entries in the data subsets via one common 
reflectivity-to-data model matrix. The optimization can then be carried out within 
the multiple measurement vector (MMV) formalism described in [16, 9, 23, 22]. 

The MMV formalism is used for solving matrix-matrix equations for an unknown 
matrix variable whose columns share the same support but have possibly different 
nonzero values. We show in this paper how to reduce the synthetic aperture imaging 
problem to an MMV format. The columns of the unknown matrix are associated 
with the discretized spatial reflectivities for different directions and frequencies. The 
solution of the MMV problem can be obtained with a matrix (2,l)-norm minimization 
where one seeks to minimize the norm of the vector formed by the ^2 norms of the 
rows of the unknown reflectivity matrix. The solutions obtained this way preserve 
the common support of the columns of the unknown matrix. 

This paper is organized as follows. We begin in section 2 with the formulation 
of the imaging problem. We derive the data model, describe the complexity of the 
inverse problem, and motivate our imaging approach. The foundation of this approach 
is in section 3, where we show how to reduce the imaging problem to an MMV format. 
The imaging algorithm is described in section 4 and its performance is assessed with 
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numerical simulations in section 5. The presentation in sections 2-5 uses the so-called 
start stop approximation, which neglects the motion of the receive-transmit platform 
over the duration of the fast time data recording window. This is for simplicity 
and also because the approximation holds in the X-band radar regime used in the 
numerical simulations. However, the imaging algorithm can include Doppler effects 
due to the motion of the receive-transmit platform, as explained in section 6. We end 
with a summary in section 7. 

2. Formulation of the imaging problem. The data model is described in 
section 2.1. Then, we review briefly imaging of isotropic reflectivity functions via 
migration and £i optimization in section 2.2. The formulation of the problem for 
direction and frequency dependent reflectivities is in section 2.3 

2.1. Synthetic aperture data model. In synthetic aperture imaging we usu¬ 
ally assume that the data u{s, t), depending on the slow time s and the fast time t, can 
be modeled with the single scattering approximation. For an isotropic and frequency 
independent reflectivity function p = p(y) we have 



( 2 . 1 ) 


with Fourier transform u{s^uj) given by 



( 2 . 2 ) 


Here k = uj/c\^ the wavenumber and the integral is over points y in D, the support 
of p. The model (2.2) uses the so-called start-stop approximation, where the platform 
is assumed stationary over the duration of the fast time recording window. We use 
this approximation throughout most of the paper for simplicity, and because it holds 
in the X-band radar regime considered in the numerical simulations. However, the 
results extend to other regimes, where Doppler effects may be important, as explained 
in section 6. 

The inverse problem is to invert relation (2.2) and thus estimate p(y), given 
u{sj^t) at the slow time samples Sj = (j — l)hs, for j = 1,..., Here hg is the slow 
time sample spacing. The inversion is usually done by discretizing (2.2), to obtain 
a linear system of equations for the unknown vector p of discretized reflectivities. 
The support Vt in (2.2) is not known, so the inversion is done in a bounded search 
domain y on the imaging surface, assumed flat. We call y the image window. The 
reconstruction of p in 3^ is a solution of the linear system, as we review briefly in 
section 2.2. 

The discretization of y is adjusted so that it is commensurate with the expected 
resolution of the image in range and cross-range. The range direction is the projection 
on the imaging plane of the unit vector pointing from the imaging location y ^ y 
to the platform location. The cross-range direction is orthogonal to range. It is well 
known in imaging that the range resolution is determined by the accuracy of travel 
time estimation, which in turn is determined by the temporal support of /(t). Thus, 
it is useful to have a short pulse /(t) whose support is of order 1/H, where B is the 
bandwidth. The range resolution with such pulses is of order c/B. The cross-range 
resolution is proportional to the central wavelength, which is why the emitted signals 
are typically modulated by high carrier frequencies Uq/ (27r). If L is a typical distance 
between the platform and the imaging window and A is the length of the flight path. 
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so that the platform receives waves within a cone of opening angle AjL^ the cross¬ 
range resolution is of the order XqL/A^ where = 27rcluJo is the carrier wavelength. 
We assume that Uq ^ which is usually the case in radar. 

In synthetic aperture imaging applications like SAR, the platform emits relatively 
long signals f{t) so as to carry sufficient energy to generate strong scatter returns, 
and thus high signal to noise ratios. Examples of such signals are chirps, whose fre¬ 
quency changes over time in an interval centered at the carrier frequency uJo/i‘^7r)- 
To improve the precision of travel time estimation, and therefore range resolution, 
the returns u{sjA) are compressed in time via match-filtering with the time reversed 
emitted signal [20]. Moreover, to remove the large phases and therefore avoid un¬ 
necessarily high sampling rates for the returns, the data are migrated via travel time 
delays calculated with respect to a reference point ia fhe imaging window. The 
combination of these two data pre-processing steps is called down-ramping. 

For the purposes of this paper it suffices to assume that f{t) is a linear chirp, in 
which case the Fourier transform |/(co’)P of the compressed signal has approximately 
the simple form 

I/Ml « |/M)|l[a;„-,rB,c^„+7rB]M, (2-3) 

where i^) denotes the indicator function of the frequency interval [ 001 ^ 002 ]- The 

down-ramped returns are 

j dt' u(^s,t - t' + 2T{s,yo)'j = y ^ , (2.4) 

and we let d be the vector of the samples of its Fourier transform 

d = > d{s,u}) = (2.5) 

The size of the vector d is NgN^. 

The linear relation between the unknown reflectivity vector p and the down- 
ramped data vector d follows from (2.5) and (2.2). We write it as 

Ap = d, (2.6) 


where the entries in p G are proportional to p(yg), with y^ the Q discretization 
points of the image window and with the constant of proportionality taken to be 
the area of a grid cell. The reflectivity p is mapped by the reflectivity-to-data matrix 
A G to the data d. The assumption of frequency independent reflectivity 

leads to a set of decoupled systems of equations A{uJi)p = d(cc;/) indexed by the 
frequency cj/, where the entries of the Ng x Q matrices A{uji) are 


Aj^q{uJl) 


\ 2iuJi [r{sjSq)-risjSo)] 

(47r|r(5j) -yg|)2 


(2.7) 


Here ki = uji/c, / = !,..., j = 1,..., A^, and g = 1,..., Q. 

2.2. Imaging isotropic reflectivities. Imaging of the isotropic reflectivities 
amounts to inverting the linear system (2.6). When this system is underdetermined, 
there are two frequently used choices for picking a solution: either minimize the 
Euclidian norm of p or its £i norm. The first choice gives 

p = aM, 
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where A’*' is the pseudo-inverse of A. If A is full row rank, A’*' = A*(AA*)“^. The 
inversion formula (2.8) also applies to overdetermined problems, where p is the least 
squares solution and A’*' = (A*A)“^A*, for full column rank A. The choice of the 
imaging window y and its discretization is an essential part of the imaging process 
and, depending on the objectives and available prior information, we may be able to 
control whether the system (2.6) is overdetermined or not. We explain in Appendix B 
that by discretizing y in steps commensurate with expected resolution limits we can 
make the columns of A nearly orthogonal. This means that in the overdetermined 
case A* A is close to a diagonal matrix. We also shown in Appendix B that in the 
underdeter mined case, for coarse enough sampling of the slow time s and frequency cj, 
the rows of A are nearly orthogonal, and therefore AA* is close to a diagonal matrix. 
Thus, in both cases, At is approximately A* up to multiplicative factors, and we can 
therefore image the support of p with A*d. This is in fact the migration formula 
(1.1) written in the Fourier domain, up to a geometrical factor, since the amplitude 
in (2.7) is approximately constant for platform trajectories that are shorter than the 
imaging distance and for bandwidths B cOq. 

If we know that the imaging scene consists of a few strong, localized scatterers, 
as we assume here, a better estimate of p is given by the optimization 

min||p||i such that ||Ap — d ||2 < e. (2.9) 

Here e is an error tolerance, commensurate with the noise level in the data, and || • ||i 
and 11-112 are the and the Euclidian norm, respectively. We refer to [1, 18, 4, 14, 13] 
for studies of imaging with optimization. The main result in this context is that 
when there is no noise so that e = 0, the reflectivities are recovered exactly provided 
that the inner products of the normalized columns of A are sufficiently small. An 
extension of the optimization to nonlinear data models that account for multiple 
scattering effects in y^ is considered in [5]. A resolution study of imaging with 
optimization is in [2]. 

2.3. Imaging direction and frequency dependent reflectivities. In gen¬ 
eral, backscatter reflectivities are functions of five variables: the location y G A’, the 
unit direction vector ni and the frequency uj. Hence, 

p = p{y,m,uj). (2.10) 


This means that the down-ramped data model is more complicated than assumed 
in equations (2.2) and (2.4) or, equivalently, after discretization, in (2.5)-(2.7). In 
integral form it is given by 


d{s^uj) = f(uj)u{s^uj)e 


r exp 

= fc2|/(w)p / ciyp(y,m(s,y),w)- 

Jn 


2iw[r(s,y) - 'r(s,yo)] 

(47r|r(s) -yl)^ 


( 2 . 11 ) 


where m(5, y) is the unit vector pointing from the platform location r(5) to y in the 
image window y. In discretized form we still have a linear system like (2.6), except 
that now p is a vector of QNgN^ unknowns, the discretized values of p in the image 
window y. 

Extending the inversion approaches described in the previous section to this model 
means inverting approximately the matrix A with a very large number of columns. 
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Fig. 3.1. Schematic of the geometry for one suh-aperture centered at the location r(sg) of the 
receive-transmit platform. 


We cannot expect the migration formula (1.1) to give an accurate estimate of the 
reflectivity as a function of five variables, as pointed out in the introduction. The 
optimization approach works, but it becomes impractical for the large number 
QNgN^ of unknowns. Moreover, it does not take into account the fact that the 
entries in p indexed by the slow time and frequency pairs (j, /), with j = 1,..., 
and / = !,..., refer to the same locations on the imaging grid. 

The imaging approach introduced in this paper gives an efficient way of estimating 
direction and frequency dependent reflectivity functions of strong localized scatterers 
in y. It uses an approximation of the model (2.11), motivated by the expectation 
that the backscatter reflectivity should not change dramatically from one platform 
location to the next and from one frequency to another. Instead of discretizing p 
over all five variables at once, we discretize it only with respect to the location in 
the image window 3^, for one probing direction and frequency at a time. To do so, 
we separate the data over subsets defined by carefully calibrated sub-apertures and 
sub-bands, and freeze the direction and frequency dependence of the reflectivity for 
each subset. The grand optimization is divided this way into smaller optimizations 
for Q unknowns, which are then coupled by requiring that the unknown vectors share 
the same spatial support in the imaging window y. 

3. Reduction to the Multiple Measurement Vector framework. We present 
here an analysis of how we can write the linear relation between the direction and 
frequency dependent reflectivity and the data as a linear matrix system 

AX = D, (3.1) 

where the unknown is the matrix X with Q rows. The entries in the rows correspond 
to the discretization of this reflectivity at the Q grid points in y. Each column of 
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X depends on the reflectivity at the backscattered direction defined by the center of 
a sub-aperture and the center frequency of a sub-band. The data are segmented over 
Ma sub-apertures and A //3 sub-bands and are grouped in the matrix D. The objective 
of this section is to describe the data segmentation and derive the linear system (3.1), 
which can be inverted with the MMV approach as explained in section 4. 

We begin in section 3.2 with a single sub-aperture and sub-band. We show in 
Lemma 3.1 that with proper calibration of the sub-aperture and sub-band size, the 
reflectivity-to-data matrix has a simple approximate form. Its entries have nearly 
constant amplitudes while the phases depend linearly on the slow time and frequency 
parametrizing the data subset. This simplification allows us to transform the linear 
system via coordinate rotation to a reference one, for all data subsets, as shown in 
section 3.3. The matrix A in (3.1) corresponds to the reference sub-aperture and 
sub-band, and the statement of the result is in Proposition 3.2. 

3.1. The sub-aperture and sub-band segmentation. We enumerate the 
sub-apertures by o = 1 ,... and denote by s* the slow time that corresponds to 
their center location r(s*). The choice of the sub-aperture size a is important, and 
we address it in the next section. For now it suffices to say that it is small enough so 
that we can approximate it by a line segment, as illustrated in Figure 3.1. The unit 
tangent vector along the trajectory, at the center of the sub-aperture, is denoted by 
tc, and the platform motion will be assumed uniform, at speed Vto,. The unit vector 
from the reference location in the image window to r(s*) is ma We call it the 
range vector for the a sub-aperture. The range (distance) to the imaging window is 

La = \r{sl) - yo\. (3.2) 


Each sub-aperture is parametrized by the slow time offset from s*, denoted by 


As = s - s* e 


a a ' 

w w. ■ 


(3.3) 


We do not index it by a because it belongs to the same interval for each sub-aperture. 
The discretization of As is at the slow time sample spacing and there are 


Tig 


a 

Vhg 


+ 1 


sample points, where a/(Vhg) is rounded to an integer. Similarly, we divide the 
bandwidth in J\fj 3 sub-bands of support b < centered at and let Aoo be the 
frequency offset 


A(jO — UJ — (jj'p G 


7 r 6 , nb 


(3.4) 


We sample the sub-band with points. 

The reflectivity dependence on the direction and frequency is denoted by the 
superscript pair (tr,/d), and by discretizing it with the Q points in y we obtain the 
vector of unknowns g C^. It is mapped to the data vector with entries 

given by the samples of d{s'^ -f As^uj^ -h Auj). The mapping is via the x Q 

reflectivity-to-data matrix described in Lemma 3.1. 

3.2. Reflectivity-to-data model for a single sub-aperture and sub-band. 

Here we explain how we can choose the size of the sub-apertures and frequency sub¬ 
bands so that we can simplify the reflectivity-to-data matrix. The calibration depends 
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on the size of the imaging window J’, which is quantified with two length scales 


and 


= max |(yq-yo)-ma|, 

<7=1,.--Q 

(3.5) 

rt = max |P„(y,-yo)|. 

<7=1,■■■Q 

(3.6) 


Here Pq, = / — is the projection on the cross-range plane orthogonal to 

and / is the identity matrix. The length scale gives the size of y viewed from the 
range direction liic,, and Yj- is the cross-range size. 

The first constraints on the aperture a and the cross-range size Y^ of the imaging 
window state that they are not too small, and thus imaging with adequate resolution 
can be done with the data subset. Explicitly, we ask that for all a = 1,..., A/'a? 


> 








> 


> 1. 


(3.7) 


The inequalities on the left involve two Fresnel numbers a^/ {XoL^) and {Y^Y/(XoL^)^ 
whose magnitudes define the imaging regime. If these numbers were small, we would 
be in a Fraunhofer diffraction regime, with approximately planar wavefronts on the 
scale of the sub-aperture and of the size of the imaging window. We consider a Fresnel 
diffraction regime, where these numbers are larger and we can get better resolution 
of images. The cross-range resolution is XoLa/a^ and naturally, the middle inequality 
in (3.7) says that the image window is larger than the resolution limit. In the range 
direction we suppose that 


> ^ » Ao, (3.8) 

where c/6 is the range resolution for the sub-bands, and we used that b < B ujo- 
While we would like to have a and b large so as to get good spatial resolution 
of the unknown reflectivity, we recall that p is frozen in our discretization in the 
small frequency sub-band and in the narrow cone of opening angle of the order a/Lc^, 
with axis defined by the center r(5*) of the sub-aperture and the reference point yo- 
The larger a and b are, the coarser the estimation of the direction and frequency 
dependence of p. The more rapid the variation of p with direction and frequency, the 
smaller a and b should be to represent it, at the expense of resolution. 

There is also a trade-off between resolution and the complexity of the inversion 
algorithm. By constraining a and b so that 


Y^ 




« 1 , 


(3.9) 


and 


A 

KLl 


«1, 


KLl 


« 1 , 


(3.10) 


we can simplify the mapping between the reflectivity and the data subset, as stated 
in Lemma 3.1. This simplification allows us to use the efficient MMV framework to 
solve the large optimization problem for the entire data set, by considering jointly the 
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smaller problems for the segmented data in an automatic way. The key observation 
here is that the unknown reflectivities for each data subset share the same spatial 
support. This is what the MMV formalism is designed to capture. 

The next lemma gives the form of the reflectivity-to-data matrix in the linear 
system 




(3.11) 


for the (a, P) data subset. It is an approximation of the system (2.6) restricted to the 
rows indexed by the slow times in the a—aperture and the frequencies in the 
/3—band. The expression of is derived in appendix A. 

Lemma 3.1. Under the assumptions (3.7)-(3.10), and with the pulse model (2.3), 
the matrix eonsists of n^j bloeks indexed by the frequeney offset 

Auji, for I = 1,... ,nuj. Eaeh bloek is an Ug x Q matrix with entries defined by 


exp I - 2i{kf) + Awi/c)m„ ■ 


(47rL„)2 
-2ikpV Asj- 


jAyq -i. 


Ln 


+ ikp - 


jAyq 


Ln 




(3.12) 


where k/, = oj*^/c, and Ay, = Yq - Yo- 

3.3. Multiple sub-aperture and sub-band model as an MMV system. 

It remains to show how to write equations (3.11) in the matrix form (3.1) with a 
reflectivity-to-data matrix independent of the sub-apertures and sub-bands. This is 
accomplished via a rotation, that brings all the sub-apertures to a single reference 
sub-aperture. But to do this, we need to know that each data subset has a similar 
view of the image window. Mathematically, this is expressed by the following two 
additional constraints on a and b 


max 

i<a<AfocA<q<Q 



mi) ■ Ay,I <C 1, 


(3.13) 


and 


max 

l<a<^^c,,l<h<^fp3<q<Q 





(3.14) 


The constraint (3.13) states that the imaging points remain within the range resolution 
limit b/c for all the apertures. The constraint (3.14) states that the imaging points 
remain within the cross-range resolution limits, as well. 

The derivation of the linear system (3.1) is in appendix A and the result is stated 
in the next proposition. 

Proposition 3.2. Under the same assumption as in Lemma 3.1 and in addition, 
supposing that eonditions (3.13) and (3.14) hold, we ean eombine the linear systems 
(3.11) in the matrix equation (3.1). The referenee sub-aperture and sub-band are 
indexed by a = 1 and 13 = 1. The unknown matrix X has Q rows and AfaJ^y eolumns 
indexed by {a,/3). Its entries are 


q 




exp 


2ikgmc, • Afq + ikg 


Afq • F^Afq 
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(3.15) 










where 


(3.16) 


rq 


p{yq,ma,U}l), 


rrin 


r(Sa) - Yo 

|r(s*) -Yol' 


The data matrix D has UsU^ rows andMaMjs eolumns indexed by (a, P). We organize 
the equations in bloeks indexed by the frequeney AuJi, for I = 1,... The entries 

of D are defined in terms of the down-ramped data veetors as 


Df'^\Au:i) = (3.17) 

where we reeall that 

d^^'^\Asj,Auji) = 6 /( 5 * + Asj.uj'^ + Auji), ( 3 . 18 ) 

and d{s^uj) is defined in (2.5). The refieetivity to data matrix A has bloeks indexed 
by Auji, denoted by A{Auji). Eaeh bloek is an Us x Q matrix with entries 


Aj^q{AuJi) = exp 


^.Auji ^ ^ 

2i -mi • Ayq 

c 


... . 

2iki—-—-t 


1 • - 


lAy^ 


(3.19) 


Note that the product of the reflectivity-to-data matrix A with each column of X 
can be interpreted, up to a constant multiplicative factor, as a Fourier transform with 
respect to the range offset mi • Ay and cross-range offset ti • Pi Ay in y. Equation 
(3.15) shows that the columns of X differ from each other by a linear phase factor 
in Ay, which amounts to a rotation of the coordinate system of the a sub-aperture, 
and a quadratic factor which corrects for Fresnel diffraction effects. Thus, the linear 
system (3.1) gives roughly the Fourier transform of the reflectivity p for different range 
direction views, and the imaging problem is to invert it to estimate p. 

4. Inversion algorithm. Here we describe the algorithm that estimates the 
reflectivity by inverting the linear system (3.1). By construction, the columns of the 
Q X AfaAfp unknown matrix X have the same spatial support, because they represent 
the same spatial reflectivity function. Thus, we formulate the inversion as a common 
support recovery problem for unknown matrices with relatively few nonzero rows 
[19, 6, 9, 12]. This Multiple Measurement Vector (MMV) formulation has been studied 
in [12, 6, 19] and has been used successfully for source localization with passive arrays 
of sensors in [16] and for imaging strong scattering scenes, where multiple scattering 
effects cannot be neglected, in [5]. 

In the MMV framework the support of the unknown matrix X is quantified by 
the number of nonzero rows, that is the row-wise Iq norm of X. If we define the set 

rowsupp(X) = {g = 1,..., Q s.t. ||e^X ||£2 7 ^ 0}, (4.1) 

where e^X is the g—th row of X and is the vector with entry 1 in the g—th row 
and zeros elsewhere, then the row-wise £o norm of X is the cardinality of rowsupp(X), 

Fo(X) = I rowsupp(X)|. 

To estimate X we must to solve the optimization problem 

minFo(X) s.t. AX = D, 
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(4.2) 







but this is an NP hard problem. We solve instead the convex problem 

minJ 2 p(X) s.t. AX = D, (4.3) 

which gives, under certain conditions on the model matrix A [ 12 , 5], the same solution 
as (4.2). In (4.3) J 2 ,i denotes the (2, l)-norm 

m 

J2,i(X) = ^||eJX||,„ (4.4) 

q=l 

which is the ii norm of the vector formed by the £2 norms of the rows of X. Fur¬ 
thermore, because data are noisy in practice, we replace the equality constraint in 
(4.3) by ||AX - D|| i? < e, where || • \\f is the Frobenius norm and e is a tolerance 
commensurate with the noise level of the data. 

There are different algorithms for solving (4.3) or its reformulation for noisy data. 
We use an extension of an iterative shrinkage-thresholding algorithm, called GeLMA, 
proposed in [17] for matrix-vector equations. This algorithm is very efficient for solving 
^i-minimization problems, and has the advantage that the solution does not depend 
on the regularization parameter used to promote minimal support solutions, see [17] 
for details. 


Algorithm 1 GeLMA-MMV 

Require: Set X = 0, ^ = 0, and pick the step size /i and the regularization param¬ 
eter 7 . 

repeat 

Gompute the residual S = Y> — AX 
X^X + /iA*(^ + £:) 

e^X ^ sign(||eJX||,, - MT) 9 = 1, • • •, Q 

Z ^ Z ^ jS 
until Gonvergence 


After estimating X with Algorithm 1 , we recover the discretized direction and 
frequency dependent reflectivity using equation (3.15), 


P 


(a,/3) _ y(a,l3) 
q ^^q 


exp 


2ik^mo, • Ayg - ik^ 


Ayq • PaAyg 


(4.5) 


for the imaging points fq = fo Ay^ indexed by q = 1 ,..., Q, the sub-apertures 
indexed by o = 1 ,..., A/'a and frequency sub-bands indexed by /3 = 1 ,..., A// 3 . 

5. Numerical simulations. We begin in section 5.1 with the numerical setup, 
which is in the regime of the GOTGHA Volumetric data set [3] for X-band persistent 
surveillance SAR. Then we present in sections 5.2 and 5.3 the simulation results. 

5.1. Imaging in the X-band (GOTCHA) SAR regime. The numerical 
simulations generate the data with the model (2.2), for various scattering scenes. The 
regime of parameters is that of the GOTGHA data set, where the platform trajectory 
is circular at height H = 7.3km, with radius R = 7.1km and speed V = 70m/s. 
The signal f{t) is sent every 1.05m along the trajectory, which gives a slow time 
spacing hs = 0.015s. The carrier frequency is uJo/{27r) = 9.6GHz and the bandwidth 
is B = 622MHz. The waves propagate at electromagnetic speed c = 3 • lO^m/s, so the 
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wavelength is Ao = 3.12cm. The image window y is at the ground level, below the 
center of the flight trajectory, and the distance from the platform to its center is 
L = 10.18km. It is a square, with side length Y = Y-^ of the order of 40m. The size 
of the sub-apertures is a = 42m and the width of each sub-band is 6 = 5/15. 

Given these parameters, the nominal resolution limits are 


AoT/a = 7.56m, c/5 = 7.23m. 

The image window y is discretized in uniform steps h = 2m in range and = Im in 
cross-range, and the reflectivity is modeled as piecewise constant on the imaging grid. 
The image discretization affects the quality of the reconstruction with £i optimization. 
It must be coarse enough so that uniqueness of the £i minimizer holds, and yet 
fine enough so that modeling errors due to off-grid placement of the unknown are 
controlled. We refer to [2] for a study of this trade-off. 

To illustrate the performance of the algorithm, we present in the next two sections 
results for various imaging scenes consisting of small scatterers supported on one pixel 
of the imaging grid, or over multiple adjacent pixels. The latter is for representing 
larger scatterers for which the direction dependent reflectivity can be motivated by 
Snell’s law of reflection at their surface. 

The results presented in the next sections compare the images obtained with 
reverse time migration and the algorithm proposed in this paper, hereby referred to 
as the MMV algorithm. The migration image is computed with the formula 


1(9) = 


(4^)^ 

kolfii^oWnsU^hh-^ 


ris rioj 

EE d{sj,uji)\v{sj) 

i=i i=i 


y|V 


■2iwi [T(sj,y)-T(sj,yo)] 


(5.1) 


which is a weighted version of (1.1), where the weights are chosen so as to provide 
a quantitative estimate of the unknown p. That is to say, when we substitute the 
data model in (5.1), under the assumption of an isotropic and frequency independent 
reflectivity we get that X(y) peaks at the true location of the scatterers and its value 
at the peaks equals the true reflectivity there. 

Let us verify the assumptions (3.7)-(3.10) with the GOTCHA parameters. The 
Fresnel numbers are larger than one, as stated in (3.7), 

^ = 5.55 and = 5.04. 

XqL XqL 

The size of the imaging region and the range resolution satisfy (3.8). Moreover, 


b Y^ 
ujo XqL I a 


0.0036, 


which is consistent with (3.9), and (3.10) is satisfied as well. 


a^Y^ 


a‘^Y 


0 . 022 . 


5.2. Single frequency results. We begin with imaging results at the carrier 
frequency, where we assume we know the range of the scatterers and seek to recon¬ 
struct their reflectivity as a function of cross-range and direction. The image window 
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Fig. 5.1. Estimation of an isotropic, frequency independent refieetivity as a funetion of eross- 
range, using No. = 8 conseeutive, non-overlapping apertures. The exaet refieetivity is shown with 
the full green line, the migration result with the blue line and the MMV inversion result with the 
broken line. The abscissa is cross-range in meters. 


extends over 120m in cross-range, and it is sampled in steps = Im, where we recall 
that XoLja = 7.56m. 

The first result displayed in Figure 5.1 is for an isotropic, frequency independent 
reflectivity of 11 scatterers, = 8 consecutive, non-overlapping apertures and noise¬ 
less data. We display in green the true reflectivity, in blue the reflectivity estimated 
with formula (5.1), and with broken line the result of the MMV inversion algorithm. 
In the legend we abbreviate the migration formula result with the letters KM, stand¬ 
ing for Kirchhoff Migration. The figure shows that the MMV algorithm reconstructs 
exactly the reflectivity, and that the weighted migration formula (5.1) does indeed 
give quantitative estimates of the reflectivity. However, the migration estimates de¬ 
teriorate when the reflectivity is anisotropic and frequency dependent, as illustrated 
next. 

The results displayed in Figure 5.2 are obtained with Ma = 10 consecutive, non¬ 
overlapping apertures. The reflectivity depends on two variables: the cross-range 
location and the scattering direction, parameterized by the slow time s*, for a = 
1,..., 10. In discretized form it gives a matrix T^true index corresponding 

to the pixel location in the image window, and column index corresponding to the sub¬ 
aperture. The reconstruction of this matrix is denoted by 7^. The green and broken 
lines in the top plots in the figure display the true and reconstructed reflectivity at the 
peak direction, vs. cross-range. Explicitly, for each pixel in the image i.e., each row q 
in T^true display the maximal entry. The migration image of the reflectors is 

independent of the direction and is plotted with the blue line. The results show that 
we have 6 small scatterers, which are well estimated by the MMV algorithm even for 
noisy data. The migration method identifies correctly the locations of the 6 scatterers, 
but the reflectivity value is no longer accurate because only a few sub-apertures see 
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Fig. 5.2. Estimation of the reflectivity as a function of direction and cross-range location for a 
scene with 6 scatterers. The top plots show the reflectivity as a function of cross-range (the abscissa 
in meters), for the peak directions. The left plot is for noiseless data and the right plot is for data 
contaminated with 10% additive noise. The green line is the exact peak value and the broken line 
the one obtained with MMV. The blue line is obtained with migration. The bottom plots display the 
reflectivity of each scatterer as a function of sub-aperture i.e., the slow time index a = 1,..., 10^ 
where 10 is the number of sub-apertures. The left plot is for the true reflectivity, the middle plot is 
for the noiseless reconstruction and the right plot is for the noisy reconstruction. 


each reflector, as we infer from the bottom plots described next. This also implies 
a deterioration in the cross-range resolution which is more visible in the next set of 
results in Figure 5.3. Naturally, the migration image gives no information about the 
direction dependence of the reflectivity. 

In the bottom plots in Figure 5.2 we show the value of the reflectivity of each 
scatterer as a function of direction, parameterized by the slow time s*. That is to say, 
we identify first the row indexes q in T^true ^ which we have a strong scatterer 
(see top plots) and then display those rows. The left plot is for the true reflectivity, the 
middle is for the noiseless reconstruction, and the right is for the noisy reconstruction. 
We observe that the MMV method reconstructs the direction dependent reflectivity 
exactly in the noiseless case, and very well in the noisy case. 

In Figure 5.3 we illustrate the effect of the anisotropy of the reflectivity on the 
imaging process. We display the results the same way as in in the previous figure. 
The point is to notice that while the MMV method estimates accurately the direction 
dependent reflectivity in all cases, the migration method performs poorly when the 
anisotropy is strong, meaning that each scatterer is seen only by one sub-aperture at 
a time (top plots). The resolution is not that corresponding to the actual aperture of 
10a = 420m, but that for a single sub-aperture of a = 42m. The middle and bottom 
row plots show how migration images improve when the anisotropy of the reflectivity 
is weaker and more sub-apertures see each scatterer. 

5.3. Multiple frequency results. Now we consider multiple frequency sub- 
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Fig. 5.3. Imaging results for anisotropic reflectivities, Ma = 10 sub-apertures and data con¬ 
taminated with 10% additive noise. From top to bottom we decrease the anisotropy. This can be 
seen from the right column plots which show the reflectivity of each scatterer for each sub-aperture. 
The middle column shows the reconstructed reflectivity as a function of direction with the MMV 
algorithm. 


bands and thus seek to estimate the reflectivity as a function of range, cross-range, 
direction and frequency. We have A/L sub-bands of width 6, and we sample each of 
them at = 15 frequencies. The number of sub-apertures is A/’a = 8- The imaging 
region is a square of side 40m and it is sampled in cross-range in steps = Im and in 
range in steps h = 2m. We denote, as before, by T^true matrix of discretized 

reflectivities and by 7^ the reconstructed ones. These are matrices of size Q x 
and we display them in the image window y as follows: For each pixel in the image 
window i.e., a row q in T^true display the maximum entry, the peak value 

of the reflectivity at point over directions and frequencies. Once we identify the 
location of the scatterers from these images, i.e., determine their associated rows, we 
display the entries in these rows, to illustrate the direction and frequency dependence 
of their reflectivity. These are the middle and right plots in the figures. 

We begin in Figure 5.4 with a single frequency sub-band (A/^ = 1), = 8 

consecutive, non-overlapping sub-apertures and data contaminated with 20% additive 
noise. The anisotropic reflectivity model has four scatterers, as illustrated in the top 
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Fig. 5.4. Results for a single frequency sub-band and Ma = 8 consecutive, non-overlapping sub¬ 
apertures. On the top we show the true reflectivity as a function of location (middle) and direction 
(right). On the bottom we show the reconstructed reflectivity with 20% additive noise. Left plot is 
the migration image, middle plot is the MMV image and the right plot is the directional dependence 
of the reflectivity reconstructed with the MMV algorithm. The axes in the left images are cross-range 
and range in meters. The abscissa in the right plots are sub-aperture index a = 1,... ,Moi. and the 
ordinate is the index of the scatterer (from 1 to 4). 


plots. Each scatterer is seen by a single sub-aperture. The reconstructed reflectivity 
is shown in the bottom plots. On the left we show the migration image, which is 
blurry and is unable to locate the weaker scatterers. The MMV algorithm gives an 
excellent reconstruction as shown in the middle and right plots. 

The results in Figures 5.5 and 5.6 are for A/L = 8 consecutive, non-overlapping 
frequency bands and flfa = 8 consecutive, non-overlapping sub-apertures. The differ¬ 
ence between the figures is the strength of the scatterers and their anisotropy. The 
results in Figure 5.5 show that the MMV algorithm reconstructs well the location of 
the scatterers and the direction dependence of their reflectivity. The frequency de¬ 
pendence of the weaker scatterers is not that accurate, likely because the bandwidth 
is small and all frequencies are similar to the carrier. As in Figure 5.4, the migration 
image is blurrier and does not locate the weak scatterers. Figure 5.6 shows that the 
migration image improves when all scatterers are of approximately the same strength 
and they have weaker anisotropy. 

The last illustration considers a larger scatterer with direction dependent reflectiv¬ 
ity, supported over four adjacent pixels, and a small isotropic scatterer with frequency 
dependent reflectivity. The data is contaminated with 20% additive noise. We note 
that the migration method gives the correct location of the large scatterer, but not the 
value of its reflectivity. Moreover, it gives a blurry image of the small scatterer. The 
MMV algorithm determines well the support of both scatterers, as well as accurate 
estimates of the reflectivity as a function of direction and frequency. 
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Fig. 5.5. Results with A/L = 8 frequency intervals, Ma = 8 apertures and data contaminated 
with 20% noise. On the top we show the true reflectivity as a function of location (left) and direction 
(middle), and frequency (right). In the middle row we show the reconstructed reflectivity with the 
MMV algorithm. The plot in the bottom row is the migration image. The axes in the left images 
are cross-range and range in meters. The abscissa in the middle and right plots are the sub-aperture 
and sub-band index and the ordinate is the index of the scatterer (from 1 to A). 


6. Doppler effects. All the results up to now use the start-stop approxima¬ 
tion of the data model, which neglects the motion of the platform over the fast time 
recording window. Here we extend them to regimes where Doppler effects are im¬ 
portant. We begin in section 6.1 with the derivation of the generalized data model 
that includes Doppler effects, and an assessment of the validity of the start-stop ap¬ 
proximation. Then we explain in section 6.2 how to incorporate these effects in our 
imaging algorithm. 


6.1. Data model with Doppler effects. For simplicity we first derive the data 
model for an isotropic reflectivity p = p(y). Then we extend it in the obvious way to 
direction and frequency dependent reflectivities in a sub-aperture indexed by a and 
sub-band indexed by /d, with reflectivity 
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Fig. 5.6. Results with A/L = 8 frequency intervals, Mot. = 8 apertures and data contaminated 
with 20% noise. On the top we show the true reflectivity as a function of location (left), direction 
(middle), and frequency (right). In the middle row we show the reconstructed reflectivity with the 
MMV algorithm. The plot in the bottom row is the migration image. The axes in the left images 
are cross-range and range in meters. The abscissa in the middle and right plots are the sub-aperture 
and sub-band index and the ordinate is the index of the scatterer. 


The scattered wave u{s,t) recorded at the transmit-receive platform is given by 


i{s,t) = - f rfydp f dti f dt 2 f”{t 2 )G{ti-t 2 ,r{s + t 2 ),y)G{t-ti,y,r{s + t)), 

Jq Jo Jo 


= -\ f dyp{y) 

Jn 




(47r)2|r(s + t 2 (t)) -y||r(s + i) -y| 


( 6 . 1 ) 


where t 2 {t) is the solution of the equation 


|r(5 + t 2 ) - y\ ^ ^ _ |r(s + ^) -y| 


( 6 . 2 ) 
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Fig. 5.7. Results with A/L = 8 frequency intervals, Mot. = 8 apertures and data contaminated 
with 20% noise. On the top we show the true reflectivity as a function of location (left), direction 
(middle), and frequency (right). The reflectivity is supported in 5 pixels, four of them are adjacent 
and represent an extended scatterer, which is frequency independent. The other pixel supports a 
small scatterer with frequency dependent reflectivity. In the middle row we show the reconstructed 
reflectivity with the MMV algorithm. The plot in the bottom row is the migration image. The axes 
in the left images are cross-range and range in meters. The abscissa in the middle and right plots 
are the sub-aperture and sub-band index and the ordinate is the index of the pixel in the support of 
the scatterer. 


and we used the expression of the Green’s function of the wave equation 


G{t, f, y) 


(5[t- |r-y|/c] 
47r|r-y| 


and the single scattering approximation. The expression (6.1) is simply the spherical 
wave emitted from r(s + ^ 2 ), over the duration ^2 of the pulse, scattered isotropically 
at y, and then recorded at v(t + s). Up to the single scattering approximation, this 
is an exact formula. Expanding with respect to t the arguments in (6.1) and (6.2) we 
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obtain 


u{s,t) = -\ [ dyp{y) 

Jq 

f" (t(l- 7 (s,y) + 0(-^-^))-2T(s,y))/(l + 7 (s,y) + 0(-^-^)) , (6.3) 
where we introduced the Doppler factor 7 defined by 


(47r|r(s) -y\y{l + 0{Vt/L)) 
VVt'- 


7(s,y) 


r'{s) 

C 


m(s,y), 


m(s,y) 


r(s) - y 

|r(s) - r{y)\' 


(6.4) 


We assume that the platform is moving at constant speed V along a trajectory with 
unit tangent denoted by t(s), and with radius of curvature R assumed comparable to 
the range L. Thus, 


7(s,y) = O(^) « 1, 

because the platform speed is typically much smaller than c, the wave speed, and 
we can neglect the residual in (6.3) which is even smaller than 7 , because over the 
duration of the fast time window the platform travels a small distance compared with 
the radius of curvature Vt R ^ L. We have thus the data model 


1 f 

u{s,t )^—^ / dyp{y)- 
^ jQ 


t{l - 2j{s, y)) - 2 r(s, y) (l - j{s, y)) 


(47r|r(s) -y|)2 


(6.5) 


which includes first order Doppler effects. 

The start-stop approximation is valid when the Doppler factor in the argument 
of /" in (6.5) is negligible. Although 7 is small, /" oscillates at the carrier frequency 
ujo which is large and, depending on the scale of the fast time t, the Doppler factor 
may play a role. Recall that t is limited by the slow time spacing hg. In practice the 
duration of the fast time window may be much smaller than although it must be 
large enough so that the platform can receive the echoes delayed by the travel time, 
2T(s,y). Explicitly, 


t = 0{L/c) + 0{l/B), 


where L/c is the scale of the travel time and 1/B is the scale of the duration of the 
signal. 

We conclude that the start stop approximation holds when 


= 0(^7) + 0(^7) « 1. 

In the GOTCHA regime, considered in the numerical simulations in section 5, we 
have 


uJqLV 
c c 


0.469, 


B c 


2.3- 10“^ 


so y) is slightly less than one. We may include it in the data model, but it 

amounts to a constant additive phase that has no effect in imaging. To see this, let 
us take the Fourier transform with respect to t in (6.5) 


S(s, uj) 


/ (^yp(y)/[w (1 + 27 ( 5 , y)) 
Jn 
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exp [ 2 i:^(l+ 7 (s,y))r(s,y)] 
(47r|r(s) -y|)2 


( 6 . 6 ) 















and expand the arguments over the slow time s and imaging point y. We use the 
approximation 


r'(s) « V t(s*) - n(s*) 


FAs 


i? J 


(6.7) 


where As is the slow time offset from the center s'^ of the aperture, and t(s'^) is the 
unit tangent to the trajectory of the platform at the center point. The second term in 
(6.7) accounts for the curved platform trajectory, with unit vector n(s*) orthogonal to 
t, in the plane defined by t and the center of curvature, and R the radius of curvature. 
We also have 


|m(s,y) -m(s*,yo)| =o(b^ 



and 


^^oT{s,y) = ujoT{s*,yo) + 0{koVAs) + 0{koAy). 

Substituting in ( 6 . 6 ) and using the parameters of the GOTCHA regime, we see that, 
^^7(3, y)r(s, y) « ujo'lis*, y^T(s*, y^), 

SO indeed, the Doppler effect amounts to a constant phase term. 

6.2. Imaging algorithm with Doppler effects. The model of the down- 
ramped data with the Doppler correction follows from (6.6), 

d{s, Lo) = /[w(l + 27 (s, yo))] w(s, u) exp [ - 2iui{l + 7 ( 3 , yo))T(s, y^)] 

« fc^/[u;(l + 27 (s,yo))] [ dy /[w(l + 27 ( 5 ,y))]p(y)x 
Jn 

exp [ 2 zu)(l+ 7 (s,y))r(s,y) - 2 iu;(l + 7 (s,yo))T(s,yo)] 

( 47 r|r(s) -yl)^ 

We are interested in direction and frequency dependent refiectivities, so to use formula 
(6.8), we consider next the o—th sub-aperture and the ^—th sub-band, where we can 
replace p by p^^’^^(y). The data is denoted by where As = s — s'^ and 

Auj = uj — Up. The goal of the section is to include Doppler effects in the statements 
of Lemma 3.1 and Proposition 3.2, which are the basis of our imaging algorithm. 

We begin with the observation that 

W 7 (s, y)r(s, y) =-^ ■ (r(s) - y) = uj'y{s, yo)T(s, y^)-Ay, (6.9) 

c c c c 

where Ay = y — yo, and f'(s) is given by (6.7), and assume henceforth that 

VY-^ b 

—« 1 . ( 6 . 10 ) 

c Lq, Uq 

This is consistent with our previous assumptions because and V and 

allows us to approximate the Doppler factor in the argument of the Fourier transform 
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of the signal in (6.8) by its value at the reference point. Then, using equation (2.3) 
and noting also that 



r / a \ / YJ- \ 1 


r /6 m 



o 

II 

l + o(-) 


we can simplify the amplitude factor in (6.8) as 


fc2/[w(l + 27(s,yo))]/[w(l + 27(s,y))] _ 

(47r|r(s)-y|)2 ^ (47rL„)2 ’ 


( 6 . 11 ) 


and obtain 




kllficOoW 






_2,(fc, + Al)Q2i±M.Ay,+ 


(47rL„)2 ^ 

2i(w^ + Aw) [t(s* + As, fo + Afg) - t(s* + As, y^)] 


( 6 . 12 ) 


Here we have used that k = -h A/c, with center wavenumber = uj^^jc and offset 
Ak = Aujjc. 

The difference between the travel times in the phase in (6.12) is approximated 
in the proof of Lemma 3.1 in appendix A. It remains to expand the first term in the 
phase, which is due to the Doppler factor. We use (6.7) and obtain 


{k /3 + Ak) -Ay, = kp— • Ay,-— 

C Cl it 


with negligible residual under the assumption 

V a 

-— <C 1. 

cRc/b 


Ha ■ Ayj + Ak^ta ■ Ay,+ 

/V a Hg- Ay, \ 
\cR cjh 1 


O 


(6.13) 


Recall that c/6 is the range resolution, and although we want cjb, the inequality 

(6.13) is easily satisfied because a R and V c. 

The generalization of the result in Lemma 3.1 is as follows. We have the linear 
system of equations 




(6.14) 


where the reflectivity vector with entries is mapped to the data vector 

with entries AuJi) by the reflectivity-to-data matrix The 

entries of are given by 


^;“’«(Awo = 


-2i 


fcol/K)l' 

(47rL„)2 
.kpVAs 


exp 


I - 2i{ki3 + Acoi/c) : 


A - ^ - A - 

xAyq - “a ■ 


1 

Ay, + —ta ■ Ay,J 


■ iki 


Ay, ■ Pa Ay, 


(6.15) 
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The difference between this reflectivity-to-data matrix and the one given by (3.12) in 
Lemma 3.1 comes from the V dependent terms in the square brackets in the phase, 
due to the Doppler effect. 

We extend next the statement of Proposition 3.2. We proceed as in appendix A, 
and show that the matrix-matrix equation (3.1), AX = D, still applies, with the same 
definition (3.17) of the data matrix D, 

and with the unknown matrix 


^(«,/3) = p{»,0) exp I _ 2ikp ■ Ay, + m„ ■ Ay 


■ ik. 


Ay, ■ Pa Ay, 


|. (6.16) 


This is under the assumptions that 


Vb 

max - 

l<a<JVoc,l<q<Q C C 


[ta - ti] ■ Ay, 


« 1 , 


V a 


l<a<Afa,l<q<Q C XqR 


max 


|(na - ill) • Ay,| < 1, 


(6.17) 

(6.18) 


which are similar to (3.13)-(3.14), and easier to satisfy for smaller V. The expression 
of the entries of the reflectivity-to-data matrix is a simple modification of that in 
equation (3.19), 


Aj^q{AuJi) = exp 


- (mi ■ Ay, + ^ti ■ Ay,) 

- 2iki (ti • Pi Ay, - L . Ay,) 


(6.19) 


Thus, the problem can be solved with the MMV approach, as described in section 
4. The Doppler correction has two effects: It gives an extra rotation in the cross-range 
direction of the imaging window (the first phase term in (6.16), involving V), and two 
extra phase factors (involving V inside the parentheses) in the reflectivity-to-data 
matrix A in (6.19). 

7. Summary. We have introduced and analyzed from first principles a synthetic 
aperture imaging approach for reconstructing direction and frequency dependent re¬ 
flectivities of localized scatterers. It is based on two main ideas: The first one is 
to segment the data over subsets defined by carefully calibrated sub-apertures and 
frequency sub-bands, and formulate the reflectivity reconstruction for each subset 
as an £i optimization problem. The direction and frequency dependence of the re¬ 
constructed reflectivity is frozen for each data subset but varies from one subset to 
another. The second idea is to fuse the sub-aperture and sub-band optimizations by 
seeking simultaneously from data subsets those reconstructions of the reflectivity that 
share the same spatial support in the image window. This is done with the multi¬ 
ple measurement vector (MMV) formalism, which leads to a matrix £i optimization 
problem. The main result of this paper is showing that synthetic aperture imag¬ 
ing of direction and frequency dependent reflectivities can be formulated and solved 
efficiently as an MMV problem. 

Data segmentation is a natural idea that has been used before for synthetic aper¬ 
ture imaging of frequency dependent reflectivities [21, 11]. Here we use it for es¬ 
timating the direction dependence of the reflectivity, as well. We analyze how the 
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size of the sub-apertures and frequency sub-bands in the data segmentation affects 
the resolution of the reconstructions as well as the computational complexity of the 
inversion. There is a trade-off in resolution in this approach: On one hand we want 
to have large sub-apertures and frequency sub-bands to get good spatial, range and 
cross-range, resolution of the reconstructed reflectivity. But on the other hand we also 
want to have small sub-apertures and frequency sub-bands to resolve well the direction 
and frequency dependence of the reflectivity. Small sub-apertures are also desirable 
so as to get images efficiently using Fourier transforms. The MMV formalism that we 
have introduced in this paper, and the associated algorithm for its implementation, 
deal well with these issues, as indicated by the numerical simulations. 

Nearly all synthetic aperture imaging is done with reverse time migration algo¬ 
rithms, without regard to whether the reflectivities that are to be imaged are direction 
dependent or not. If the reflectivities are isotropic, then the spatial resolution of the 
reconstruction improves as the aperture increases. But this is not the case with 
direction dependent reflectivities as only part of the synthetic aperture will sense 
reflectivities from particular locations. This means that segmenting the data over 
sub-apertures is natural. The MMV-based imaging algorithm introduced in this pa¬ 
per handles automatically signals received by sub-apertures that are coming from 
directional reflectivities located in the image window. 
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Appendix A. Derivation of the reflectivity to data model. Here we show 
that the expression of Aj^q{uJi) in (2.7) can be approximated by given in 

Lemma 3.1, for cj/ = -1- Acj/ and Sj = s* -h Asj. Eor simplicity of notation we drop 

the indexes j and I of the frequency and slow time. 

It is easy to see from (2.3) and the assumptions ujo'^h and Lq, a > IG that 


(47r|r(5) — y|)^ 


(A.l) 


ioT k = uj/c and ko = uJq Ic. It remains to show the phase approximation 



, Ay.P^Ay 

- T - 


(A.2) 


where cj = -h Acj lies in the frequency sub-band of width 6, s = 5* -h As is in the 
sub-aperture of size a and y = yo + Ay is in y. 

We begin by expanding the travel time in Ay, 


$ = 2w [T(s,y) - r(s,yo)] 
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with small residual 








by assumption (3.10) and < a, inferred from (3.7). Here we used the expression 
of the gradient 


Vy|r(s) -y| 


r(g) - y 
|f(s) - y| 


= -m(s,y), 


the Hessian 


' Vj;|f(s) -y| = 


l?(s) - y| L' 


J-m(s,y)m^(s,y) 


and 


bl,<7=i 

Next, we expand in Auj = oj — uj'^ and obtain 

ki3 


3ma ■ Ay 

|r(s) - y|2 


[|Ay|^ - (m„ • Ay)^]. 


$ = -2{kis + Afc)m(s, fo) ■ Ay + 


|f(s) - fo 


Ay ■ [/ - m(s, yo)m^(s, y^)] Ay + £’ 2 , 


where Ak = Aujjc and 






< 1. 


The last estimate is by assumption (3.9). Finally, we expand in As = 5 — s*, and 
recalling the notation in section 3.1, we get 


^ = —2{k^ + Ak)ma, • Ay — 2/c; 


HAs- 

' T ^ 

-L/fy 


»Ay+ 


, Ay-P„Ay , ^ 

fc/3- f - \-£- 


(A.3) 


The residual is the sum of four terms 


£ — £^2 T £*3 T £*4 T £^51 


with S 2 given above. The term £3 comes from the quadratic part of the expansion of 
fc/ 3 m(s,yo) • Ay, 


£3 - kpiVAsf 


-ng ■ PgAy 

RLa 


tg ■ [ifiginj + mg(mg)^] Ay 

VL^ 


(tg ■ PgAy) (tg ■ mg) 1 

LI \ 


Here ~ denotes order of magnitude, and the primes denote derivative with respect to 
s. The unit vector rig is normal to tg, in the plane defined by tg and the center of 
curvature of the trajectory of the platform. It enters the definition 


t' 


Vng 

R 


(A.4) 


where R ^ is the radius of curvature. Moreover 

















We conclude that 


^3 = 0( 


KLl 



«1, 


where the inequality is by assumption (3.10). 
The term £ 4 ^ in the residual is 




Acj VAs - 


.Ay = 0( 


r b aYg \ 

\ UJn ^nLrv ' 




«1, 


by assumption (3.9), and the last term £^ comes from the expansion of the quadratic 
term in Ay in the expression of We estimate it as 


^5 


V XoLl 


) 


< 1 , 


where we used assumption (3.10). The statement of Lemma 3.1 follows from (A.l) 
and (A.3). □ 

Proposition 3.2 follows easily from the expression (3.1) of and assumptions 

(3.13) and (3.14). Writing the linear system (3.11) component-wise we get 


q=l 


exp 


AuJi ^ ^ 1 

— 2i -• Yq — 2ikfsVAsj- 






(Awi), 


with given in (3.15) and defined in (3.17). The result (3.19) follows 

from this equation and assumptions (3.13) and (3.14). □ 

Appendix B. Inner products for rows and columns of the reflectivity- 
to-data matrix. 

Here we analyze the relation between the discretization of the imaging window 
y and the linear independence of the columns of the reflectivity to data matrix. 
This is done by computing inner products of of normalized rows and columns of the 
reflectivity-to-data matrix. If the column inner products multiplied by the number 
of elements in the support of the reflectivities are below a threshold then the MMV 
algorithm will give an exact reconstruction, in the noiseless case [5] . 

We consider the restriction to a data subset, defined by a sub-aperture and fre¬ 
quency sub-band satisfying the assumptions in section 3. Thus, we work with matrices 
A(a,/5)^ but to simplify notation we drop the indexes (<a,/3). 

Let us denote by the q—th column of matrix A and calculate the inner product 


(ug,a^/) — ^A A^ 

\ / q,q‘ 


ris 

EE ^j' 5<7(A^0 (Acij/). 


i=i 1=1 


Using Lemma 3.1 we get 


= exp 1 ^ - 2 iki3ma, ■ {yq' - jq) + 


jfc/3(Ay,/P„Ayg/ - AygP„Ayg^ 




EE 

j=i 1=1 


exp 


2iAuJi ^ ^ . 2ikRVAsj^ ^ 

-■ (y,, - fq) -A-J t« ■ Pa 


(yg' - Yq) 
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where we normalized the columns by their Euclidian norm. The sums can be approxi¬ 
mated by integrals over the frequency band and aperture, as long as they are sampled 
at intervals h^j and hg satisfying 

huj |ma ■ (yg - yg')| Vhs \ta ■Fa{yq> - yq)\ 

b c/b ’a XoL/a ' 

We obtain after taking absolute values that 


sine -hIq 
\c 


iyq' - yg))sinc(^ta • Pa(yg' - y,)) 


(B.l) 


This is small for q q' when we sample the imaging window y in steps that are larger 
than the resolution limits c/b in range and XoL/a in cross-range. 

A similar calculation can be done for the rows of A, denoted by We have 




and using Lemma 3.1 we get 




1 ^ 


q=l 


exp 


2i(cj// -UJi)ma- Afq 


q=l 


2ik^V{sj' - Sj)ta • IPaAyg 


Furthermore, for discretizations of the imaging window in steps h in range and in 
cross-range, satisfying 

\uJi'-LOi\ h V\sj-Sj^\ 

b c/b ' a \oLa/a ’ 

we can approximate the sum over q by an integral over the imaging window and obtain 

Hj'.n 


{uJl' UJi)Yq,\ . / koV{Sj' 


\||a(,.,,,)|n|aov,,)||/ \ c ) \ 


La 


This result shows that the inner product of the rows is small when the frequency is 
sampled in steps larger than Y^/c and the slow time is sampled in steps larger than 
{l/V)/{XoY^^/La,). 
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